Coupled electron— heat transport in nonuniform thin film semiconductor structures 
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A theory of transverse electron transport coupled with heat transfer in semiconductor thin films is 
developed conceptually modeling structures of modern electronics. The transverse currents generate 
Joule heat with positive feedback through thermally activated conductivity. This can lead to insta- 
bility known as thermal runaway, or hot spot, or reversible thermal breakdown. A theory here is 
based on the optimum fluctuation method modified to describe saddle stationary points determining 
the rate of such instabilities and conditions under which they evolve. Depending on the material 
and system parameters, the instabilities appear in a manner of phase transitions, similar to either 
nucleation or spinodal decomposition. 

PACS numbers: 72.60.+g, 72.80.Ng, 64.60.Q-, 73.50.Fq 
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I. INTRODUCTION 

Various treatments of electronic transport in disor- 
dered systems typically concentrate on systems at a given 
fixed temperature. However, observations (see references 
below) often point at the coupled electron-heat trans- 
port where local fluctuations in electric current generate 
temperature fluctuations. When the latter have positive 
feedback, as e. g. in the case of thermally activated con- 
ductivity, an instability arises leading to the current fil- 
amentation. 'Weak spots' corresponding to suitable dis- 
order configurations promote such instabilities. While 
this mechanism has long been known qualitatively^ its 
quantitative understanding remains insufficient leaving 
open questions about the role of material and structure 
parameters, and effects of static vs. thermodynamic fluc- 
tuations. 

This work attempts a theory of coupled electron-heat 
transport concentrating on a rather representative case 
of transverse conduction through thin-film structures. 
A model structure consists of an active (heat generat- 
ing) conducting layer between two electrically inactive 
insulating layers representing encapsulation always found 
with electronic devices. This structure is depicted in 
Fig. [1] The active layer can be a single or multi- 
layered semiconductor sandwiched between thin metal 
electrodes. The electric potential along each of the elec- 
trodes is constant; the potential difference V between 
them is maintained by an external power source. 

The disorder is introduced through the activated 
transversal electric conduction with random Gaussian ac- 
tivation barriers varying in the lateral (along the film) di- 
rections. The role of insulating layers is that they affect 
the temperature distribution and make the entire model 
more realistic. For simplicity, we assume one of them to- 
tally insulating while another one having a finite thermal 
conductivity. Also, for simplicity, thermal conductivities 
and specific heats of the active and insulating layers are 
assumed the same. 

The analysis below is aimed at finding the probability 
of local temperature fluctuations and their radii associ- 
ated with locally increased current density vs. the sys- 
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FIG. 1: Sketch of the system with nonuniform power genera- 
tion and current flow (fat arrows) . Top insulating layer shown 
in gray. The bottom layer (dark) represents a strong thermal 
insulator. 



tern dimensions, material parameters, and ambient tem- 
perature. It is based on the premise of localized rare 
lateral fluctuations that do not overlap. These local- 
ized entities are similar to other types of localized states 
in disordered systems, for which theoretical description 
known as the optimum fluctuation method (OFM) has 
been developed long ago. OFM was originally created 
to describe electronic states in band tails of disordered 
semiconductors j^r— it was applied later to localized sound 
excitations in glasses^, resonance electronic states in 
disordered metals^^ fluctuation tail states in magnetic 
semiconductors, — random lasing in disordered dielec- 
tric films jii, local fluctuations in thermal expansion of 
glasses^ and nucleation in disordered medial 

The essence of OFM is in the optimization of config- 
urational probability (or entropy) of fluctuations under 
the additional condition that the dynamical characteris- 
tic of a fluctuation satisfies the appropriate differential 
equation (Schrodinger equation for electronic state, elas- 
tic wave equation for sound excitations, electromagnetic 
wave equation for optic modes, etc.). This is achieved 
through the variational approach, in which the dynamical 
characteristic is kept fixed (yet arbitrary) in the course of 
optimization of the configurational entropy, after which it 
is optimized to additionally minimize that entropy. The 
details of OFM vary between different systems. Here de- 
veloped OFM is tailored to describe the temperature flue- 
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tuations coupled with the electric current, so that the dy- 
namical characteristic (temperature) of fluctuations sat- 
isfies the heat transfer equation. 

The analysis below shows that hot spot instabilities 
evolve in a manner of phase transformations, either by 
nucleation or similar to spinodal decomposition affecting 
the entire area. The nucleation scenario of such instabil- 
ities in uniform systems was established earlier based on 
general phenomenological analysis^ 

This paper is limited to a general theoretical anal- 
ysis; possible applications of the coupled electron-heat 
transport will be presented in more appropriate journals. 
We refer to a recent monograph [l4| for many practically 
important cases. The relevant observations are found 
with bipolar transistors^—, other metal-insulator- 
semiconductor structures^ - — , nanoscale transistors^, 
graphene transistors,—, and thin-film photovoltaics.— ~— 
In these applications, the phenomenon under consider- 
ation was labeled as thermal runaway, or hot spot, or 
(reversible) thermal breakdown. It can be detrimental 
to the corresponding device operations leading to their 
irreversible degradation in hot spots via local shunting, 
burning, or melting; hence, significance for device relia- 
bility. 

The paper is organized as follows. Sec. |TT] introduces 
the basic equations describing the coupled electron-heat 
transport in a non-uniform system. To better explain the 
essence of OFM and subsequent results, two toy models 
are considered in Sec. IIV1 Sec. |Vj presents a modifica- 
tion of OFM describing saddle points through which the 
system evolves into thermally non-uniform state. The 
OFM functional is optimized in Sec. IVII through direct 
variational procedure. The steady state rate of hot spot 
nucleation is estimated in Sec. I VIII Finally, Sec. IVIIII 
presents general discussion and conclusions. 



II. COUPLED ELECTRON AND HEAT 
TRANSPORT IN A DISORDERED SYSTEM 



The Joule power density is given by 



E 



P = P exp(-P//cT), Po = ^oexp^--j. (1) 

Here £ = V/ho is the electric field strength where ho is 
the distance between the electrodes (see Fig. [T]). <7 is 
the pre-exponential of conductivity, 

<r = a cxp[- (E + E)/kT] 

with E being the average activation energy, k is the 
Boltzmann's constant, and T is the local temperature. 
The random part of activation energy, E has zero aver- 
age, (E) = and a finite dispersion (E 2 ) = B. It is 
characterized by the correlation function 



Here the radius vector r lies in the film plane, S(r) is 
the two-dimensional delta function implying zero correla- 
tion radius disorder. The minimum area s is determined 
by the physical nature of fluctuations. For example, its 
characteristic linear scale ao ~ s 1 / 2 (likely in sub-micron 
range) can be given by the screening radius or the grain 
size, or other length, below which the system parameters 
do not vary significantly, s is introduced to give B the 
dimensionality of the square of energy and the meaning 
of the dispersion of random energies E. 

Local elements of the system interact through heat 
transfer described by the standard equation 



X V 2 T 



P(r) = 



(3) 



where x is the thermal conductivity and the Laplacian 
V 2 is three dimensional, and x is coordinate independent. 
The power generation density is a sum of average and 
random contributions, 



P= (P) 



P (i) 



(P) = Pn 



b \ exp 



E 
kf 



where 



P (1) = P exp 



E 
kf 



(P). 



(4) 



Eq. ([3|) assumes the steady state heat transfer. The 
assumption of stationary states is common to all known 
cases of OFM. The problem under consideration, how- 
ever, is different with respect to the notion of stationary 
fluctuations. Since the instability evolves in a fashion of 
phase transitions, the stationary solutions of Eq. ([3]) can 
only describe saddle points in the parameter space. The 
temperature fluctuation ST becomes time dependent in 
the proximity of each of such point, described by 



-CST/t = xV 2 T + P(r) 



(5) 



in the relaxation time approximation, where C is the spe- 
cific heat. The fluctuation decay will correspond to pos- 
itive, while fluctuation growth (instability) to negative 
values of r; this criterion is used in Sec. IVII below. 



III. LINEAR APPROXIMATION: 
NO-BREAKDOWN STEADY STATE REGIME 

For completeness, consider briefly a trivial situation 
where the disorder B and temperature fluctuations ST 
are small allowing the linearization 



P = Po 



1 + 



E(r) 
kT 



(6) 



where Tq is the average temperature. Substituting this 
into Eq. §3§ and setting 



(P(r)P(r')) =BsS(r-r'). 



(2) 



ST(r, z) = 4>{r) exp(z/zo), zq = const 
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for the radial (r) and transversal (z) coordinates yields 

1 



V#- -j0 = «(r). 

Here V 2 is the two-dimensional Laplacian, 

1 _ PoS = X kTg P Q E(r) 



(7) 



4P* 



xkT 



and Zq must be determined from the boundary condi- 
tions. The solution to Eq. ([7]) has the form 

0(r) = (-1/4) J d 2 ru(v')H^\i\v - r'|/r ) (8) 

where Hq is the Hankel function. 

The quantity in Eq. ((8j) represents a sum of large num- 
ber of random contributions and, according to the central 
limit theorem, is a random quantity itself with the Gaus- 
sian probability distribution. Its dispersion (</> 2 ) is given 
by 




4x 2 (fcT ) 2 4 x Ek ' 

Here we have taken into account Eq. ([2]) and the valued 
of the integral / [Hq(x)] 2 xcIx — 2. 

To avoid unnecessary discussions of boundary condi- 
tions zo and E are left as two parameters. Neglecting the 
temperature change through the active film (exp(/io / z ) ~ 
1), the net result is that the temperature fluctuations are 
characterized by the radii of ro and the Gaussian distri- 
bution, 



ST 2 

p(ST) oc exp ( --^2 



with S T 2 = ^g^. (10) 
4 X Ek 



The important point is that the above linear approx- 
imation does not account for positive feedback of tem- 
perature fluctuations on transversal conduction and thus 
the disorder remains fixed and temperature independent. 
While this restriction eliminates the possibility of ther- 
mal breakdown (which is the main topic here), the results 
of this section can still be applicable to the case of very 
small currents and fluctuations used e. g. in thermogra- 
phy diagnostics 



IV. TOY MODELS 

Because the regular OFM below is mathematically 
cumbersome, it is illustrated here with simplified (toy) 
models. One of them concentrates on the case when there 
is no positive feedback on conductivity by local heating. 
Another one deals with a homogeneous system and con- 
centrates on the positive feedback. 



A. Conductive filaments through an insulating film 

Consider a two phase structure where transversal cur- 
rent flows through conductive filaments in an insulating 
host of thickness ho sandwiched between two equipoten- 
tial electrodes. The structure is characterized by the av- 
erage transversal conductivity a due to filaments of av- 
erage concentration n per area. Local fluctuations Sn 
in their concentration result in the corresponding con- 
ductivity fluctuations Sa = aSn/n. Since the filaments 
generate Joule heat, they create fluctuations ST in tem- 
perature; the tail of probabilistic distribution of ST is 
found below. 

Consider a cylinder shaped region of radius a perpen- 
dicular to the electrodes where the characteristic fluc- 
tuation in filament concentration is Sri. The Gaussian 
probability of such a fluctuation is estimated as 



exp 



(Sn) 2 a 2 



exp 



-no 2 



exp(-S) 



S can be optimized with respect to a after Sa is expressed 
via ST and a.. 

The heat flux through the cylinder base and side sur- 
faces is estimated as x[(5T / h )a 2 + (ST / a)h a\. Equating 
it to the fluctuation of power V 2 a 2 Sa/h inside the cylin- 
der yields the temperature fluctuation 



ST = 



V 2 Sa 



X 



(12) 



Expressing Sa from Eq. (1121) and substituting it into Eq. 
(HH yields 



S = 



hi 

a z 



(13) 



Following the OFM approach, we optimize the expo- 
nent S with respect to the fluctuation radius a, i. e. 
dS/da — 0, which gives a = ho. Substituting a = ho 
back into Eq. (fT3| yields the optimum exponent of prob- 
ability, 



Sopt — 



ST 
ST 



where STq = 



V 2 a 
xhoVrl 



(14) 



again to the accuracy of numerical multipliers. 

The preexponential is roughly estimated by dividing 
the entire area into elemental domains of area /iq each 
and noticing that exp(— S op t) describes the probability 
of a desired fluctuation with temperature excess ST in 
a given domain. Therefore, the concentration of such 
fluctuations is estimated as h^ 2 exp(— S op t). 

Two features should be noted. First, OFM concen- 
trates on the exponent of probability, largely neglect- 
ing the pre-exponential factors (although they can be 
estimated as well). Secondly, it optimizes that expo- 
nent in order to find the most likely disorder configu- 
ration providing the desired fluctuation characteristic of 
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interest. Its applicability is limited to the region of non- 
overlapping fluctuations. 

A possible application of this toy model can might be 
a system of multiple shunting metal chains formed in 
dielectric or solid electrolyte films considered for non- 
volatile memory; see e.g. Refs. [36| and references therein. 



B. Homogeneous Alms 

Consider, in the linear approximation, a relatively 
small temperature fluctuation ST in a cylinder region of 
radius a, setting 



1 _ 1 ST 

T~y ~T$- 



(15) 



Neglecting (for simplicity) heat transfer through the 
cylinder bases and using 

r _ fSTE\ 
Sa = ae * P [-kT*)> 



(16) 



Eq. (|T2j) reduces to the form 



V 2 aa 2 (5TE\ 
SI = — exp 



X 



h 2 ~* \kT 2 J ' 



For a system in equilibrium, the probability of tem- 
perature fluctuation ST in volume SV — ira 2 ho is given 
by the expression^ exp[-C^SV(ST) 2 /kT 2 ] where 
is the specific heat per volume. Expressing a 2 from 
Eq. (1161) gives the equilibrium distribution function 
J(5T) oc exp[-S(<JT)] with 



S(ST) = - 



2kT 2 V 2 a 



ST 1 



exp 



STE\ 



(17) 



It follows from Eq. (fTTj) that the equilibrium dis- 
tribution is a minimum at ST C — kT 2 /3E where the 
product ST 3 enp(—STE/kT 2 ) is a maximum. This can 
be interpreted as a barrier in the system free energy at 
ST = 5T C : the probability of fluctuations first exponen- 
tially decreases as ST grows below ST and then decreases 
when ST exceeds ST C . Such a behavior is obviously sim- 
ilar to that known in nucleation phenomena^— (where 
the barrier is a function of the nuclear radius) and small 
polaron collapse^ (where the barrier is a function of di- 
lation). The instability point corresponds to a relatively 
very small temperature increase ST C = (kTo / '3E)Tq <C Tq 
in systems with high enough activation energies, say, 
ST C < O.OITq - 3 centigrade. 

Based on that analogy, the exponent of probability of 
the thermal breakdown is given by S(ST C ), that is, to the 
accuracy of numerical multipliers, 



S(ST C ) 



k 2 T%C^hl x 
E 3 V 2 a 




(18) 
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FIG. 2: Effective barrier for nucleation of hot spots corre- 
sponding to the numerical value a — EkTo = 10. Arrows 
show a pathway of hot spot nucleation. 



Note that the probability exponent optimization here 
results not in a minimum, but rather a maximum; it may 
turn into a saddle point in a parameter space of higher 
dimensionality as will be explicitly shown next. Another 
conclusion is that a positive feedback alone makes the 
instability possible regardless of the degree of disorder in 
the system. 



V. OPTIMUM FLUCTUATION METHOD 

The subtlety of the optimum fluctuation method is in 
how it treats the disorder induced distribution of tem- 
perature T(r) (or wave function for the standard case of 
energy spectra in systems with random potential energy). 
Namely, T(r) is considered a smooth 'optimum' function 
approximating the temperature distribution for the most 
likely disorder configuration responsible for any desired 
temperature fluctuation. It remains arbitrary (yet fixed) 
in the course of the analysis and is determined later by 
the condition of the maximum of the probability. Such 
optimization benefits from the known property of varia- 
tional techniques that any inaccuracy in the trial function 
translates into a higher order inaccuracy in the corre- 
sponding functional. 

In what follows we take into account only exponen- 
tially strong activation factor ignoring all possible pre- 
cxponcntials found with temperature dependent conduc- 
tivity in semiconductors. This simplification simultane- 
ously determines the accuracy of our analysis where all 
the pre-exponential factors are replaced with their aver- 
ages. In particular, this analysis is limited to the case of 
strong enough fluctuations beyond the linear approxima- 
tion for P[6T(r)]. 
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A. OFM equations 

The heat transport equation ([3]) can be treated as an 
extremum of the functional 



F= d 



(VT) 2 - P(r) 



(19) 



where the pre-exponential factor (E + E)/T 2 [generated 
by variation of P in Eq. ([T])] is approximated by its 
average, 



C = X ((E + E)/kT 2 ). 
The latter functional can be presented as 



F 



d 3 r 



(P) 



Z 



(20) 



(21) 



where T depends on coordinates and random variable Z 
is defined by 



Z = J d 3 rP (1) (r) 



(22) 



OFM suggests that the dispersion of random variable 
Z can be found as 



D = (Z A 



d 3 rd 3 r'{P {1 \v)P {1 \v')) (23) 



where the average in the integrand is evaluated under the 
condition of a fixed (yet arbitrary) function T(r). The 
integral in Eq. (I22p contains a large number of random 
contributions. Therefore, according to the central limit 
theorem, Z is described by Gaussian statistics, i. e. its 
probabilistic distribution 



j(Z) cc exp[-S(Z)], S=—. 



(24) 



The maximum probability fluctuation corresponds a 
stationary point of S(Z) under the additional condition 
of Eq. (|21[) . Finding such a conditional extremum is 
tantamount to finding an unconditional extremum of a 
functional 



Z 2 
2D 



XF 



(25) 



where A is the undetermined Lagrange multiplier. A is 
then found from the additional condition of a certain pre- 
determined maximum temperature in the the optimum 
fluctuation region. 

The functional $ must be optimized with respect to 
the disorder configuration E(r) and the field T(r). Be- 
cause the former appears only with the integral Z, the 
optimization can be more conveniently conducted with 
respect to Z and T(r). The corresponding equations are 



Z 
D 



A = 



(26) 



and 



Z 2 SD 



2D 2 6T 



A£V 2 T 



k 2 T 3 



exp 



(E 2 ) 
2k 2 T 2 



= 0. (27) 



Here we have taken into account a known property 3 -^ 

(exp(P/fcT)) = cxp[((P/fcT) 2 )/2] 

for a Gaussian random variable E/kT. Using Gaussian 
statistics in combination with the concept of thermally 
activated current assumes the inequality 



E (E 2 ) 
— » - — -. 

kT k 2 T 2 



(28) 



Allowing the opposite inequality would lead to the phys- 
ically unacceptable feature that the typical fluctuation 
current exponentially decreases with temperature. 

Substituting Eq. ^ into Eqs. (JUJ) and (J27|) yields 
the equations determining the optimum fluctuation tem- 
perature field T(r) and its corresponding probability ex- 
ponent, 



XDSD _ , 



f E (E 2 )\ f (E 2 ) 
[kT 2 ' ~ fc2 T 3 ) ex P {2k 2 ~T 2 ~ 



= 0, 



S = 



DX 2 



(29) 



(30) 



To evaluate 5D/ST that is the variational derivative 
of the integrand in Eq. (123[) we use again the property 
of averaging of a Gaussian random variable E(r). The 
integrand in Eq. (f2"3"|) becomes 



Pq exp 



f (E 2 ) - 




_(kT) 2 _ 


J d 3 r' jexp 



(E(r)E(r')) 



k 2 T(r)T{r') 



- 1 



For the case of delta correlated disorder in Eq. <j2j> , the 
latter expression can be approximated as 



Pgs/io exp 



2B 

Jkfy 



(31) 



Substituting the result of differentiation [together with 
Eq. J2HJ)] into Eq. (|2"T|) leads to a closed form single equa- 
tion for the optimum fluctuation T(r). That equation is 
not very useful practically because of its rather complex 
form . The problem becomes easier when presented in 
the form of functional subject to direct optimization with 
respect to T(r). That functional is given by 



J = 



where 

F= |(VT) 2 -P exp 



d 3 rF[T(r)} 



B 



— XPqV exp 



2P 



(32) 



(33) 



2(kT) 2 \ u ^L(fcT) 2 

Note that, to the accuracy of the factor of —A, the third 
term in the functional J [corresponding to the third term 
in Eq. (|33[) ] is twice the probability exponent S. 
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B. OFM saddle points 

While optimization of functional J remains to be im- 
plemented, the nature of its stationary points can be 
determined already here. Assuming a trial function 
T = T(r/a) and changing variable r — > r/a, J can be 
presented in the form 



J 



J i + a 2 J 2 



proach to nucleation (Zeldovich' theory; see e. g. Chap- 
ter XII in Ref. |32| ) and in agreement with the quali- 
tative analysis in Sec. IIV B[ that minimum determines 
the nucleation barrier and rate. This approach will be 
implemented in Sec. IVIII below upon determining the 
parameters of OFM solutions 5T(r). 



VI. DIRECT VARIATIONAL PROCEDURE 



where J\ and J 2 do not depend on a. Treating a 2 as 
a variational parameter, leads to the conclusion that 
d 2 J/d(a 2 ) 2 = at the stationary points where J 2 = 0. 
Hence, they represent inflection points rather than min- 
ima. In a higher dimension parameter space including 
the temperature fluctuation amplitude, these points can 
only be saddles. 

The saddle point solutions require a different interpre- 
tation of OFM results. From the physical standpoint, 
some (but not all) of their related configurations should 
appear with certainty, i. e. with S — 0, since they are 
not steady state, and thus are to be passed inevitably 
sooner or later. From that perspective, they are simi- 
lar to the barriers of classical nucleation theory^ir— or 
small radius acoustic polaron formation^ For example, 
the OFM saddle points in the surface J (a, T) can phys- 
ically describe critical radii a(T) separating the regions 
of spontaneous decay from that of spontaneous growth 
of fluctuations. This similarity to the nucleation theory 
will be made explicit in Sec. I VII 

Note that the fact of probability exponent S vanishing 
at the OFM saddle points, does not compromise OFM as 
long as the corresponding fluctuations remain strongly 
localized and do not overlap. The latter conditions do 
not necessarily invoke S > 1 (unlike the conclusion of 
Sec. IIVI where all the fluctuations simultaneously coex- 
ist), since the saddle point events are not steady state 
taking place at different time instances. 

Consider the configurational probability exponent S 
in a certain proximity of a saddle point S = 0. We de- 
note STq(t) the temperature distribution in the optimum 
fluctuation corresponding to S = 0. If the optimum fluc- 
tuation ST(r) is different from STq(t), one can extend 

S = Jd 3 r (0s^ [ST(r) - M>(r)] 2 , (34) 

where the integrand is positive. The equilibrium distri- 
bution function of such fluctuations is given by 



/(5T) = / exp 



C(») 
2fcT 2 



d 3 r[ST(r)} 2 - S \ (35) 



Here f is the preexponential factor and we have taken 
into account the expression for the probability of equi- 
librium temperature fluctuation ST in volume SV men- 
tioned in Sec. HVBl 

It is seen from Eq. (|35[) that / is a minimum at some 
ST different from STq. Following the Fokker-Planck ap- 



A. Trial function and functional 

Here we implement a direct variational procedure of 
optimization of the functional J using the simplest trial 
function 



ST 



1 



r 

ah 



1 



when ST > 



(36) 



that is zero outside of the domain r < dh, z < h. Here r 
and z are the radial and transversal (across the film) co- 
ordinates. 9 and d are the two variational parameters, de- 
fined as being dimensionless to make the resulting equa- 
tions more compact. In particular, is the amplitude 
excess temperature in fluctuation measured in the units 
of the average temperature To, and a has the meaning of 
the fluctuation radius measured in the units of structure 
thickness h. 

Note that integration over the transversal (z) coordi- 
nate extends over the entire structure thickness (h) for 
the first term in Eq. (1331) . while the second and third 
terms must be integrated only over the active layer thick- 
ness (ho <C h) where the power is generated. Also, we 
note that the constraint t = at z = h correctly reflects 
the boundary condition of a constant temperature at the 
interface (see Fig. [IJ. Furthermore, we assume fluctua- 
tion to be relatively small, allowing the linearization in 
Eq. (H5J). 

Substituting Eq. (l36l) and carrying out the integration 
reduces J to the form 



12a 2 J 
where 



(d 2 



-f3d 2 <Z>(x)-\f3 f3d 2 <Z> x— (37) 



1 / s exp(x) — x — 1 ,„„ N 

x = a6 and $(x) = \ . (38) 

x l 

Here we have introduced the parameters defined as 



B 



E 



kT (fcT ) 2 ' 



2ihh P a 2 
p = — 5 exp 



a =2 



B 



(kT 



B 



2(kT ) 



, P = Po^exp 



2a, 



35 



2(fcT ) 



The inequality in Eq. (|28| limits them to a > 1. In 
integrating over z in Eq. (I37|) . we have assumed a prac- 
tically important case when the semiconductor layer is 
very thin, a6ho/h <C 1, and calculations are simpler. 
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Because eventually we consider 9 = x/a an in- 
dependent given variable, the optimization conditions 
d J I da 2 = and d J/ dx = must be used to solve for a 2 
and A. In agreement with the conclusion of Sec. [V] the 
stationary points found from the optimization are saddle 
points. This is seen from the sign of the determinant 



d 2 j d 2 j 

(da 2 ) 2 (d6) 2 



d 2 J 



l 2 



{da 2 )d9 



< 



identifying the stationary points as saddles^ 

B. Regional approximations 

Consider the results of optimization of the functional 
J for three complimentary regions. 

1. Weak fluctuations, x <C 1 

Assuming x <C 1 reduces $(x) in Eq. ([57]) to $(x) ps 
1/2 + x/6 + x 2 /24, which significantly simplifies the op- 
timization. This leads to the physically unacceptable so- 
lution with a 2 = -32/(8 + /3) <0. 

2. Moderate fluctuations, x ~ 1 




10 12 14 



a0 



FIG. 3: Phase diagram for a thin film structure with transver- 
sal current vs. power density (parameter /3) and local tem- 
perature increase (parameter a9). Region to the left of the 
line aO = 4 represents the stable phase where local temper- 
ature fluctuations decay making thermal breakdown impossi- 
ble. The gray colored region below the line of solution of Eq. 
(|43p . represents metastable state corresponding to the saddle 
points, through which thermal breakdown nucleates locally. 
The solid curve in that region is a solution of Eq. (|44|) ; it 
corresponds to the most likely nucleation events, for which 
S = in Eq. (139[) . The region above the line of solution of 
Eq. (|43p represents the globally unstable state of the system. 



It is straightforward to verify that the interpolation 
$(x) = 1/2 + x 2 /6 holds to the accuracy of several per- 
cent for intermediate x < 4. Using that interpolation, 
the optimization of J results in the physically inconsis- 
tent solution as well, a 2 = -[12 + 16(ev6>) 2 ]/(6 + 3/3). 



3. Strong fluctuations, i>l 

Acceptable solutions with a 2 > exist in the case of 
a9 1 (and yet a9ho/h <C 1) where one can approxi- 
mate $(x) = exp(x 2 )/x 2 . This yields 



A = 



d 2 



(a9) 4 ~ /3 exp{a9) 
W exp(2a<9) ' 
4(afl) 3 
2{a9) 4 -/3exp(a0)' 

6>[(a6>) 4 - Pey^{a9)] 2 exp(-2a6l) 



where 



and 



So 



2(ad) 4 - P exp(a6 



4er 2 ) 2 cx P [-2g/(fcr ) s 

288P 2 w/i 



9 c i < 9 < 9 C 2, 



(39) 



(40) 



(41) 



(42) 



with t c \ and t C 2 being the two solutions of the transcen- 
dental equation 



2(a9) 4 - /3exp(a6 



0. 



The condition 



(a9) 4 -(3exp(a9) = 



(43) 



(44) 



describes the points where S = and thus nucleation of 
hot spots takes place, according to the discussion in Sec. 
IV Bl These points all fall within the domain of physically 
acceptable solutions in Eq. (|42|) . Also, it follows from 
comparison of Eqs. (|43|) and (jH} that the radii of the 
corresponding stationary fluctuation states remain finite 
as required by OFM. 

Because (aO) 4 exp(— ad) is a maximum at a9 = 4, Eq. 
(l44l) has solutions when 



4.7 



(45) 



where [e] stands for the base of natural logarithms. Close 
to that threshold value, the dependence t(/3) takes the 
form 



aB a9 = 4 + ^/3 C - (i when / 8 C - J 9<1. (46) 

Another branch of a# with the minus sign before the 
square root is ignored as belonging to the moderate fluc- 
tuation regime. 
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Alternatively, one gets from Eq. (|44|) . 

aJaa9 = ln(l/j8)»l when /3 < /3 C . (47) 

This behavior corresponding to the far right part of the 
solid curve in Fig. [3] describes the low power regime. 

C. Phase diagram 

The complementary region to the left of the line a6 = 4 
in Fig. [3] was characterized by the physically unaccept- 
able solutions with a 2 < (see Sec. IVIB II and IVI B 2j) . 
Here, we argue that that region represents the state 
where the system remains stable with respect to ther- 
mal fluctuations. A proof is achieved by including in the 
above analysis the term —C5T/t from Eq. ([5|) describ- 
ing the temporal behavior of fluctuation. It is straight- 
forward to see that the unacceptable negative o 2 turn 
positive when r > 0, i. e. the corresponding fluctuations 
decay. 

Alternatively, for the region above the curve ft — 
2(a9) 4 exp(— x), adding the term with negative relax- 
ation time r < allows for positive a 2 . Therefore, the 
states in that region are globally unstable, i. e. they 
evolve into highly conductive high temperature states 
without any barrier. This is qualitatively similar to the 
phase transition scenario of spinodal decomposition,— 
which is not described in the OFM framework. 

Note the triple point O at (j) = j3 c , a9 = 4) in Fig. [3] 
where all three phases coexist. It is straightforward to 
show that fluctuations 59 become increasingly strong in 
its proximity where 



Sof3 2 a 



4) 2 (S9f 



(48) 



and \59\ = \9 — 9o\ <C 9q. That property is similar 
as well to that of the standard phase transition phase 
equilibria^ 



D. Approximation of classical nucleation theory 

The approximation of classical nucleation theory im- 
plies a narrow boundary region between the two phases 
and its related concept of surface energy. It can be at- 
tempted in the current framework by choosing a trial 
function 



ST 



1 when r < a, 
(a + d — r) / d when 
when r > a + d 



a < r < a + d, (49) 



with rf«(i. As a result, the gradient term in Eq. (|3"3")l 
is determined by the contribution from a narrow layer of 
width d analogous to nucleus interfacial energy in func- 
tional J of Eq. (|3"2"j) . The procedure of optimization 
becomes even simpler than that based on the trial func- 
tion of Eq. (|3"o| . Omitting the details, the result is that 



the functional J has no stationary points when d <C a. 
Hence, the approximation of interfacial energy does not 
apply to the case under consideration; the function in Eq. 
P6p remans more adequate. 



VII. STEADY STATE TRANSITION RATE 

Consider the probability of thermal breakdown at a 
given power density Pq described in terms of the parame- 
ter p < f3 c . Using ST(r, z) from Eq. (|36[) and expressions 
for a and S from Eq. (|39[) . the equilibrium distribution 
function becomes 



f(9) = f exp 



irC^h 2 h 9 2 d 2 (a9) 



3k 



S{a9) 



(50) 



S(a9) is a maximum, S = 0, at the line shown in Fig. 
[3] and increases towards the boundary a.9 = 4. However, 
given realistic parameters (see Sec. IVIII|) that increase is 
not nearly as significant as the increase of the first term 
in the exponent in Eq. (|50l) . As a result, f(9) has a sharp 
minimum at a9 ~ A. 

Following the known approach of nucleation theory^ 2 - 
(mentioned in Sec. IV Bl above) consider a stationary 
Fokker-Planck equation 



df 

j = -B-qq +Af = const 



(51) 



for the 'kinetic' temperature distribution function f{9). 
Here j is the flux in the temperature fluctuation (9) 
space, D is the diffusion coefficient in that space; A 
is connected with D by a relationship which follows 
from the fact that j = for the equilibrium distribu- 
tion / = /. Using the latter enables one to present 
the flux as_j = -Bf(d/d6)(f //), and, hence, f/J = 
—s J d9 1 Bf ' + const. Finally, applying the boundary con- 
ditions / — > when t — > oo and f — f when 9 = 0, 
yields 



d9 
W 



(52) 



The integral is determined by a narrow proximity of the 
minimum of / that gives the exponent of the transition 
rate. 

To roughly evaluate the preexponcntial factor (without 
any knowledge of D) one can divide the entire area into 
a set of cells of characteristic linear size of the optimum 
fluctuation ah. Then the preexponential must be of the 
order of the rate of temperature variations n/{ah) 2 in a 
cell where k is the thermal diffusivity. This yields the 
steady state nucleation rate (cm _2 s _1 ), 



3 



16k 



exp 



nC^h 2 h 9 2 a 2 (A) 



3A- 



5(4) 



(53) 



where a 2 (4) = a 2 (a9 = 4) and 5(4) = S(a9 = 4) are 
given in Eq. (I39[) The power density enters this result 
through the parameter j3 in Eq. (|4ip . 
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This result becomes more explicit for the case of low 
enough power when /3exp(4) < 4 in Eq. (|59")l and the 
absolute value of the exponent in Eq. (f53f is estimated 
as 



5w 8 



C^h 2 h Q 
a 2 k 



-7-10 



_ 6 (£T 2 ) 2 exp[-2£/(fcT ) 2 ] 



otti^sP^ 



(54) 



This is similar to the exponent in Eq. (1181) emphasizing 
the important role of specific heat and rapidly decreasing 
with the power density. However it has a distinct feature 
of a lower boundary beyond which it cannot be further 
reduced even for very high power densities. It should be 
remembered however that high enough power densities 
are conducive to a different type of instability similar to 
the spinodal decomposition transformations as reflected 
in Fig. H 



VIII. DISCUSSION AND CONCLUSIONS 



g(«0 



5T r 



ST 



FIG. 4: Probability g(8T) of hot spots vs. their excess tem- 
perature ST. The Gaussian tail at low ST is described in Sec. 
IIVBI The critical overheat ST C corresponds to the condition 
aQ = 4 illustrated in Fig. [3] The high temperature peak 
at STh is determined by the processes of saturation of acti- 
vated conduction and inter-spot interactions as explained in 
Sec. IVIIIBI its width is due to disorder effects. 



A. Numerical estimates 

Assuming the typical semiconductor values^ one gets 
X ~ 1 W/cm-grad and E/T ~ 10 — 100 for activation 
energies E ~ 1 eV and T - 100 - 500 °K. This yields 
£ ~ 1 - 100 W/cm-grad 2 , a ~ 10 - 100. 

For geometrical parameters, it is natural to assume 
ho ~ 1 /im, s r~j 1 fim 2 , and h ~ 10~ 4 — 10 _1 cm. The cur- 
rent density in the range from 1 /xA/cm 2 to 1 A/cm 2 and 
electric fields £ ~ 10 3 — 10 5 V/cm are used in many device 
operations. The corresponding power densities are in the 
range from 1 mW/cm 3 to 10 5 W/cm 3 . The fluctuation 
strengths exponent exp[— 2B / (kT ) 2 ] can be evaluated as 
~ 0.001 — 1 based on the observations of transversal cur- 
rents through nonuniform Schottky barriers and thin film 
photovoltaics4i Finally, we use the thermodynamic pa- 
rameters C ~ 0.1 — 1 J/sm 3 -°C and n ~ 0.1 — 1 cm 2 /s. 
With the above parameters, the preexponential factor in 
Eq. dS3J) is estimated as ~ 10 5 - 10 13 cm" 2 s _1 . Given 
that preexponential, the exponent in Eqs. (|53l) and (|54[) 
can be then within the range of experimentally impor- 
tant nucleation rates only for micron or sub-micron thin 
devices. Assuming greater thickness, say, h ^ 1 mm 
makes the thermodynamic term proportional to C large 
enough to practically rule out the possibility of thermal 
breakdown mechanism under consideration. 

However, semiconductor devices of modern electronics 
are often 10-100 nm thick (unless intended thermal sinks 
are used), and for them the thermodynamic fluctuation 
term in the exponent is not too large. For such structures, 
the second term in the nucleation rate exponents can 
be not terminally large for powers in the range P > 
100 W/cm 3 . Overall, this makes the above considered 
mechanism realistic for structures in submicron region. 

Finally, the minimum power density corresponding to 
the critical value of j3 in Eq. (|4"5"|) , above which the nucle- 
ation mechanism turns into that of global instability, can 



be estimated as Po > 10 11 W/cm 3 . This range of power 
density is above practically all types of modern semicon- 
ductor devices, except maybe some cases of power elec- 
tronics. 



B. Discussion 

The above consideration is limited to a basic instability 
triggered by Joule heat in combination with activated 
conduction. The instability is predicted to start under 
insignificant local overheats of several degrees. However, 
this analysis does not address the final parameters to 
which the instability can grow. 

The 'stabilized' temperature excess STh in the devel- 
oped filament (beyond the present theory framework) can 
be rather substantial. As pointed in Ref. lj, it can be- 
long in the temperature range where the activated con- 
duction saturates. That high temperature limit should 
not be mixed with the above predicted transition tem- 
perature excess, ST C w AkT 2 /E <C STh (corresponding 
to aO ~ 4), starting from which the instability evolves. 
This is illustrated in Fig. |4] 

Furthermore, it is conceivable that the steady state 
high temperature local overheat STh cannot be deter- 
mined by any extension of the present theory limited to 
nonintcracting hot spots, even if activated conduction is 
allowed to saturate. The concentration of steady state 
hot spots at STh can significantly depend on their inter- 
action. Indeed, the present theory predicts (Sec. IVIIj) 
that even at arbitrarily however low rates, the above de- 
scribed instabilities will keep developing (maybe beyond 
the practically significant time intervals) to take over the 
entire structure area. This contradictory prediction is not 
unique of the system under consideration. It is known in 
the theory of phase transition where the nucleation stage 
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is limited by various inter-nucleus interactions, such as 
competition for material, elastic stresses, etc. Similarly 
limiting interactions here will include competition of hot 
spots for the electric current, thermal fields by other fil- 
aments, etc. This analogy leads to the prediction of the 
growth and ripening stages of thermal breakdown kinet- 
ics, similar to that of the standard phase transitions^ 
a theory of such later stages of hot spot transformation 
remains to be developed. 

While not related to structural transformations, the 
predicted local temperature increase can accelerate such 
transformations leading to permanent failures in the form 
of conducting pathways. Therefore, this mechanism can 
serve as a precursor to permanent structural failures. 
From that point of view, the above results on low temper- 
ature thermal breakdowns ST C <C T point at high sensi- 
tivity of the fatal failure probability to the activation en- 
ergy of conductivity and thermodynamic variables, par- 
ticularly, specific heat, thickness, and thermal insulation. 

The role of inactive (thermally insulating) layers ex- 
ponentially reducing the thermal breakdown rates is due 
to the filament diameter increase with its length. As a 
result the thermal gradient in radial direction decreases 
suppressing the instability rate. This is consistent with 
the known practical solutions using substantial heat sinks 
attached to with submicron electronic devices in order to 
minimize their failure rates. 

A more theoretical comment is in order regarding the 
relevance of the above OFM modification aimed at 'non- 
traditional' saddle type of stationary points. The un- 
derlying motivation was to relate localized temperature 
fluctuations with other known localized states in disor- 
dered systems. However the same basic equations as 
derived in Sec. |V] could be obtained in the framework 
of instanton approach suitable for theoretical descrip- 
tion of nucleation42r— That approach would start with 
the time dependent heat transfer equation leading to 
the variational problem for the exponent of probability 
exp[— R(T, t)} where t is time and R is related to the func- 
tional in Eq. {TJJ]), R cx J F[T(t)]dt. F remains a ran- 
dom functional to be additionally optimized to maximize 
the probability. That reduces the conditional variational 
problem for R to that of unconditional extremum in Eq. 
(I2"5j) yielding final expressions of OFM in Eq. J2H]). 

The above theory has the following limitations. (1) 
The assumption of fixed voltage V across the film im- 
plying that the current / through the filament must be 
small enough, IR s h -C V where R s h is the sheet resis- 
tance of the conductive electrodes. (2) Simplification of 
uniform thermal conductivity may have noticeable quan- 
titative ramifications, yet can hardly change the qualita- 
tive predictions. (3) The approximation of ^-correlated 
disorder, according to which the transversal conductiv- 
ity must fluctuate across the distances smaller than the 
filament radius. The opposite regime of strongly cor- 
related disorder can be readily described by the above 



results reduced to the case of homogeneous structures, in 
which then consider Pq or a as a random quantity varying 
over distances greater than the filament radius. (4) The 
optimum fluctuation method per se with accuracy lim- 
ited to the probability exponent. (5) Inaccuracy of the 
direct variational procedure with a simplistic trial func- 
tion remains unknown. Based on many similar examples, 
one can expect the results to be semi-quantitatively cor- 
rect. (6) Limitation of small temperature fluctuations 
a6ho/h <C 1, remains self-consistent as long as it is con- 
sistent with the final results for 9 as it takes place in the 
above. 



C. Conclusions 

The following was shown. 

(i) Thin film semiconductor structures with activated 
transversal conduction are unstable with respect to re- 
versible thermal breakdowns in the form of hot spots and 
their related current filaments. 

(ii) The instabilities evolve in a manner of phase transi- 
tions by either nucleation (at not too high power densi- 
ties) or absolute instability similar to spinodal decompo- 
sition (above certain critical power density). 

(iii) The optimum fluctuation method can be modified 
to describe the saddle points, through which such tran- 
sitions occur. 

(iv) The instabilities start with finite local temperature 
fluctuations that are smaller than the average tempera- 
ture To by the factor of kTo/E with E being the average 
activation energy of electric conduction. The initial fluc- 
tuation radii are by the same factor smaller than the 
structure thickness. 

(v) The stable, metastable, and unstable phases of a ther- 
mally uniform system form a diagram (in variables power 
density - temperature) similar to the standard phase dia- 
grams of phase equilibria, in particular, with fluctuations 
diverging towards the triple point. 

(vi) The steady state nucleation rate of hot spots ex- 
ponentially depends on the material parameters, system 
geometry, and disorder strength. 

The author hopes that this consideration can form a 
theoretical basis to analyze system failures in various 
structures of modern thin film devices; specific examples 
will be presented elsewhere. 
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